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1 Introduction 

For the study of functional organization and reorgani¬ 
zation of the human cortex by means of electromag¬ 
netic source imaging, changes in the location and spa¬ 
tial extent of neural sources is of interest. We present 
a probabilistic approach based on the maximum like¬ 
lihood (ML) theory for modeling extended sources 
using anatomical and physiological knowledge to re¬ 
duce the ambiguity of the ill-posed biomagnetic in¬ 
verse problem. The number of sources, their location 
and extension is defined by a set of parameters which, 
due to measurement noise, is associated with a certain 
probability distribution defining the likelihood func¬ 
tion (LF). The parameter set which maximizes the LF 
is determined by using the Metropolis algorithm. 

To investigate the influence of different levels of un¬ 
correlated and correlated noise on the outcome of the 
proposed method, simulations with an activity located 
on both sides of an artificial sulcus were performed. 
The applicability of the proposed method is pre¬ 
sented as an application on experimental somatosen¬ 
sory evoked magnetic fields. Our tests demonstrate 
the applicability of the presented method in locating 
extended sources and quantifying their extension. 

2 Methods 

2.1 Mathematical basis 

Let b = {&i,..., b m } be the set of measurements 
representing the m values of the recorded magnetic 
field B at m sensors at a specific time point. In a sta¬ 
tistical framework these values can be viewed as the 
realisation of m random variables according to a 
certain probability distribution pi determined by the 
noise level and noise type at each sensor. The current 
distribution j representing the cortical activity under¬ 
lying the measured magnetic field B can be described 
as a physical system characterized by a set of param¬ 
eters 0q — {0oi? • • •, 0ok}- Thus, solving the ill-posed 
inverse problem corresponds to the estimation of the 
parameter set 0. 


The maximum likelihood method (MLM) associates 
each set of parameters 0 with a certain probability for 
the realisation b. For statistically independent mea¬ 
surements yi this value is given by the expression: 

m 

m = y[p i {y i \e). (i) 
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The function C(0) represents the likelihood or like¬ 
lihood function (LF). The maximum likelihood es¬ 
timator (MLE) 0 for the parameter 0 is the value 
of 0 which maximizes the LF. In practice it is 
more convenient to consider the logarithm of the LF, 
L(0) — ln{£(0)} called log-likelihood function (log- 
LF). Since the logarithm is a monotonic function, the 
log-LF is maximized by the same value 0. 

As many measurements are required to appropriately 
describe the covariances of dependent variables Bi , 
we are assuming a m-dimensional normal distribu¬ 
tion of independent random variables. We further as¬ 
sume that the noise distribution at each sensor has the 
same standard deviation a. Introducing this ansatz 
into Eq. 1 the log-LF takes the form: 

m 

m = -^2 

-ylna-yln(27r). (2) 

L i denotes the i -th row of the leadfield-function L 
and j(0) the current density distribution representing 
the cortical activity due to the parameter set 0. Since 
the last two terms in Eq. 2 are independent of the pa¬ 
rameter set 0, the maximization of L{0) of the LF cor¬ 
responds to the minimization of the quadratic form: 

m 

= j(0)» 2 - (3) 
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2.2 Anatomical constraints 

Anatomical knowledge of the neural origin of the 
measured activity makes it possible to adopt con¬ 
straints in order to reduce the ambiguity of the ill- 





Figure 1: Artificial sulcus from both views with ex¬ 
actly reconstructed activity. The white triangle repre¬ 
sents the center of the source (R — 10 mm). Its ridged 
edge is the result of inaccuracies in the exact distance 
reconstruction between two triangles. 


posed inverse problem. Since MEG activity is gener¬ 
ated mainly by the apical dendrites of pyramidal neu¬ 
rons predominantly oriented perpendicular to the cor¬ 
tical surface, the modeled neural activity can be con¬ 
fined to the cortical gray matter which can be modeled 
by a surface composed of contiguous triangles [1,3]. 
The activity of a single source triangle is represented 
by an equivalent current dipole oriented orthogonally 
to the triangle. Because of the lack of experimental 
data regarding the current source density across an 
active cortical area, a constant electric dipole moment 
per unit area is assumed, which implies a constant ac¬ 
tivation of the cortical tissue. Consequently, the mo¬ 
ment of a single current dipole is proportional to the 
area of the associated triangle. 

2.3 Parametrization of the activity 

The parameters 9i are chosen to unambiguously de¬ 
scribe the cortical activity and to incorporate knowl¬ 
edge of the anatomical and physiological facts under¬ 
lying the neural activity. As a first order approxima¬ 
tion we model sources as being circular, which im¬ 
plies nondisjunctive active cortical areas. The num¬ 
ber of synchronuous active sources is also assumed 
to be known. Consequently, the cortical activity is 
parametrized by two continuous parameters: i) the lo¬ 
cation of the source T and ii) the radius R of the active 
cortical area. The location of the source is defined by 
the index of the triangle which represents the center 
of the source. The triangles appertaining to the active 
source satisfy the condition that the distance between 
their centers of gravity and that of the central triangle 
of the source maximally equals R. 

2.4 Scanning the parameter space 

Since it is generally not possible to specify 9 ana¬ 
lytically, we propose to use an adapted Metropolis- 


algorithm (MA) based on a Monte Carlo Markov 
Chain (MCMC) approach [4], which scans the param¬ 
eter space according to the values of L(9), thus in¬ 
creasing the probability of finding the LM estimator 9 
within the performed sample. The algorithm is based 
on a diffusion process within the parameter space re¬ 
alized by a random walk of the source locations on the 
extracted cortical surface and concomitantly assign¬ 
ing randomly different values to the source radius. 

3 Simulations 

To develop and test the proposed method we used the 
simulated activity of an extended source located on 
both walls of an artificial sulcus (Fig. 1) consisting of 
1200 triangles with an average area of 1.1 mm 2 . Its 
shape and orientation were chosen to resemble the so¬ 
matosensory cortex, for which an experimental study 
is presented in a later section of this article. The simu¬ 
lated cortical source was defined by 195 triangles lay¬ 
ing within a radius of Ro — 10 mm around a central 
triangle Tq. Consequently, the extended source was 
parametrized by 9o — {To, Ro}. 

The distance between two triangles is determined by 
an iterative procedure which computes the shortest 
link between the centers of gravity of the correspond¬ 
ing triangles over the vertices and centers of gravity 
of all the other triangles. Thus, instead of calculating 
the euclidean 3D-distance between two triangles, the 
distance along the cortical surface is considered, bet¬ 
ter accounting for the functional organization of the 
cerebral cortex. 

The MEG analysis was performed for a whole-head 
MEG system with 151 first order axial gradiometers 
(CTF Systems Inc.) and a spherical volume conductor 
with an outer radius Ry — 90 mm. 

3.1 Influence of uncorrelated and correlated noise 

Since MEG recordings are inevitably corrupted with a 
certain level of noise which negatively influences the 
accuracy of the source reconstruction, it is important 
to assess the impact of this noise type on the MLM. 
For this purpose three levels of uncorrelated and cor¬ 
related noise were investigated (1%, 5% and 10%). 
Uncorrelated noise was simulated by adding normally 
distributed numbers with mean 0 to all 151 measure¬ 
ments bi. /The normal distribution was chosen such 
that the ratio between its variance and that of the an¬ 
alyzed signal resulted in the considered noise level. 
Uncorrelated noise stemming from spontaneous brain 


log-likelihood function 



Figure 2: Dependence of the log-likelihood function 
on the extension and the location of the source. 


activity was modeled by adding the potential of 1000 
current dipoles with uniforml random orientation sit¬ 
uated in an upper hemisphere with an eccentricity of 
70 mm around the center of the volume conductor. 
The amplitudes of the ECDs were chosen such that 
the ratio between the variance of their resulting ac¬ 
tivity at all sensors and the analyzed activity resulted 
in the considered noise level. Ten different simula¬ 
tions were performed for all three levels of both noise 
types. 

3.2 Results 

Figure 1 shows L{6) as a function of the source lo¬ 
cation and source extension for the simulated activ¬ 
ity in absence of noise for 50,000 computation steps. 
Although all regions are reached during the fit pro¬ 
cess, the MCMC algorithm scans predominantly the 

Table 1: Influence of noise expressed as the distance 
between the centers of the default and the recon¬ 
structed sources Ad. 


noise 

Ad (mm) 

type 

level 

min 

max 

mean 


i% 

0.0 

3.2 

0.7 

uncorr. 

5% 

0.0 

5.1 

2.4 


10% 

0.0 

5.1 

2.5 


1% 

0.0 

5.7 

3.7 

corr. 

5% 

3.2 

8.2 

5.0 


10% 

1.3 

10.0 

5.6 


Table 2: Influence of noise on the source expressed as 
the difference between the radii of the default and the 
reconstructed sources A R. 


noise 

A R (mm) 

type 

level 

min 

max 

mean 


i% 

0.0 

0.5 

0.1 

uncorr. 

5% 

0.0 

0.5 

0.6 


10% 

0.0 

4.5 

0.6 


1% 

0.0 

7.5 

3.0 

corr. 

5% 

0.0 

8.0 

4.0 


10% 

0.0 

9.0 

4.8 


parameter space in regions with higher LF values. 

As a measure of the accuracy of the reconstructed 
source, the distance between the central triangles of 
the modeled and the default sources Ad and the dif¬ 
ference A R between their radii were chosen. The out¬ 
come of the noise analysis is presented in Table 1 and 
Table 2. 

4 Application 

To illustrate the ML method for an experimental data 
set, we studied the somatosensory representation of 
the left-hand digits (dl-d5) of a subject on the basis of 
MEG-recordings using a whole-head system with 151 
axial first-order gradiometers (CTF Systems Inc.). 

4.1 Method 

The subject’s left hand digits received tactile stimula¬ 
tion in a random sequence either individually or syn- 
chronuously, resulting in 6 stimulus conditions (dl, 
d2, d3, d4, d5 and dl-d5). As the components mod¬ 
eled are in the 60-75 ms latency range, it is known 
a priori that activities result from SI [5, 2]. Thus, to 
reduce the ambiguity of the inverse problem, a region 
located one centimeter anterior and posterior to the 
right central sulcus of the right hemisphere (contralat¬ 
eral to the stimulated side) was used as the region of 
interest for the source reconstruction (Fig. 3). 

4.2 Results 

The analysis of the experimental data confirmed our 
hypothesis that the stimulation of larger skin areas re¬ 
sults in an extended activation of the somatosensory 
cortex (Table 3). 






















































Table 3: Formal radius of the modeled sources. 



stimulus condition 

dl 

d2 

d3 

d4 

d5 

dl-d5 

formal 
radius (mm) 

1.8 

1.0 

3.0 

2.1 

1.9 

10.0 


theoretical predictions. In the case of the simulta¬ 
neous stimulation of all five digits (dl-d5), the cor¬ 
responding source comprises the sources correspond¬ 
ing to the individual finger representations (Table 4, 
Fig. 3). Further, the application demonstrates that a 
reliable estimate of the source extension is possible 
despite the presence of correlated noise. 



Figure 3: Somatosensory representation of the sub¬ 
ject’s left-hand digits due to the analysis of the so¬ 
matosensory magnetic field evoked by synchronous 
stimulation of all five digits (condition dl-d5). 


5 Discussion 

Our simulation demonstrates that the presented max¬ 
imum likelihood method is able to assess the location 
and the size of an extended cortical source, even in the 
rather difficult case of opposing activity on the ante¬ 
rior and posterior wall of an artificial sulcus, where 
both activities tend to cancel each other when viewed 
from distant electrodes. 

The method, however, is limited by uncorrelated and 
correlated noise associated with the measurement and 
analysis of evoked magnetic fields. The simulations 
reveal that increasing noise causes a decrease in the 
accuracy of the reconstructed source as well as a sys¬ 
tematic underestimation of its extension. 

The source analysis using the MLM shows that the so¬ 
matosensory representation of the single digits have 
similar extensions. The locations of the individual 
sources are also in very good accordance with the 
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